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1. Introduction 

In a historical numerical experiment Fermi Pasta and Ulam |I] studied the following 
non-linear model of a ID crystal described by the Hamiltonian 

N (\ x 11 

H — I -mx 2 + V{ x i ~~ x i-i ~ ) 5 V(x) = -mu 2 x 2 + -kx 4 . (1) 
i=i / 2 4 

Contrary to initial expectations of Fermi et al the so-called FPU model (|I]) behaved 
in a strong contrast to the ergodic theorem of statistical mechanics, even in quite 
strongly non-linear regime when one would expect fast relaxation to statistical canonical 
equilibrium and equipartition of energy from an arbitrary initial state. Instead, FPU 
model triggered the discovery of non-linear normal modes, the so-called solitons, and 
indirectly, the whole field of (computational) non-linear dynamics. Later, FPU chain 
has been studied in the spirit of KAM theory [[J and various estimates have been given 
for the critical strength of the dimensionless non-linearity parameter kl 2 /mu 2 required 
for the motion to become globally stochastic || [| . 

Recently, FPU-type models have been studied in the context of energy transport 
and Fourier heat law in ID chains of particles |5|, ||, 0]. Quite unexpectedly, FPU-type 
models, namely Hamiltonians of the type ([!]) with non-linear inter-particle interaction 
V{x) (and no on-site potential), turned out to be anomalous heat conductors, due to 
slow power law decay of transport (velocity-velocity or current-current) time correlation 
functions. These results, namely that correlations decay universally as C(t) ~ t~ 3 / 5 , and 
correspondingly that Kubo transport coefficient diverges as n(L) ~ L 2 / 5 as a function 
of the chain length L, have later been successfully explained in terms of hydrodynamic 
arguments and mode-mode coupling theory M. 

However, a perturbative-like approach such as mode-mode coupling theory should 
break down when the inter-particle potential has more than one stable position. Indeed, 
Giardina et al [[| have performed a series of numerical experiments on the particle chain 
with an oscillating inter-particle potential V(x) = 1— cos(x), as well as with the potential 
V(x) = —x 2 /2 + x A /A, and found normal heat conduction with a clean exponential decay 
of correlations. The key mechanism, which we believe produces non-KAM-like (almost) 
hyperbolic motion of such a high dimensional Hamiltonian system, are the hyperbolic 
saddles over which pairs of particles flip from one well to another. This motivated us 
to study in this paper some fundamental dynamical properties of the simplest version 
of such a hyperbolic chain with quartic double- well inter-particle potential V(x). We 
present here some intriguing numerical results which suggest existence of a simple (but 
only approximate) Markov partition with very simple many-body symbolic dynamics of 
this particle chain. In addition, we find that the transition matrix is formally generated 
by a Hamiltonian of a certain quantum spin-1/2 chain. 

Section |2| introduces the model. In section |3] we define (approximate Markov) 
partition of phase space and obtain numerical and analytical estimates for the volumes 
of various cells. Expression for the average time between subsequent jumps among the 
cells is also obtained. In section |] we introduce Markovian description of our model 
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with the transition matrix and the flux matrix and numerically check the accuracy of 
the Markovian property in two independent ways. Section [5] represents the core of 
the present paper. We numerically calculate the transition matrix and make a formal 
correspondence between our model and a ferromagnetic Heisenberg quantum spin-1/2 
chain. Using this correspondence we derive various estimates and bounds for the decay 
rates of the time correlation functions of a class of piece-wise constant functions. 



2. Model: Inverted Fermi-Pasta-Ulam chain 



With the choice for units of mass, time and length of m, 1/u and I, respectively, and the 
new canonical coordinates = Xi — il and pi = mxi, Hamiltonian (P can be brought 
to a dimensionless form 

N 2 2 4 2 

fl = E| + ^-a V(x) = -a- + X - + ( 2 ) 

i=l 

where a = —mu 2 /kl 2 is the only dimensionless parameter left. In order for the potential 
minima to be at level (for a > 0) we have added a constant term a 2 /4 to the potential 
V(x). We can get all different orbits of the system described by the Hamilton function 
(0) just by varying one parameter a. At first it would seem that the energy E is also a 
free parameter. But by fixing all the three basic units we have also provided a natural 
unit for the energy density e = E/N. If we scale the dimensionless quantities by a 
factor £, namely (qi,pi,a) — > £ 2 Pi, £ 2 a), the Hamiltonian is simply multiplied by a 
factor £ 4 . Invariant parameter in this transformation is therefore A = e/a 2 . All different 
nonequivalent Hamiltonians can be obtained just by varying one parameter, namely A, 
or equivalently a at fixed e. Energy density can therefore be fixed without any loss 
of generality. From now on e = E/N = 1 is assumed, as well as periodic boundary 
conditions (qN+i,PN+i) = (qi,Pi)- Phase space point will be denoted by x = (q,p). 
For negative values of the parameter a the system has a form of the standard /3-FPU 
model (just in a different parametrization), but in this paper we focus on the positive 
values of the parameter a in which case we call the system Inverted FPU model (IFPU). 
IFPU model is translationally invariant, therefore besides the energy e there is also a 
second (trivial) constant of motion, namely the total momentum P = YliLiPi- m a U 
numerical calculations the total momentum has been fixed to P = 0. 

For the numerical integration of equations of motion a fourth order symplectic 
algorithm of Ref. [|H]] has been used. This integrator has been checked to be the optimal 
choice (at least for the model studied here, and at relative accuracy ~ 10~ 5 — 1CT 8 ) 
by making a careful comparison with a number of other symplectic and Runge-Kutta 
methods. 



3. Phase space partition and statistical dynamics 

If the barrier AU = a 2 /4 between the potential minima is sufficiently high the differences 
(<Zt + i(£) — qi(t)) will spend most of the time either around the left or the right minimum, 
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Figure 1. Potential V(x) for the inverted FPU model (||), 



a > 0. 



with quite infrequent jumps between the wells. As the motion around the minimum 
is approximately harmonic the more interesting part will be the transitions between 
the left and the right potential well, for each pair of neighboring particles, going over 
the hyperbolic saddle at q i+ i — g» = 0. Understanding the motion in phase space of 
a high- dimensional system is generally very difficult. Some symbolic description of an 
orbit x(t) is therefore highly desirable. In IFPU model with high potential barrier the 
choice of a partition of phase space is almost obvious: 

For each instant of time t we will be interested only in a single binary bit of 
information a; G {0, 1} for each pair of neighboring particles (i,i + 1), namely if the 
neighbors are in the right or left well, (qi+i — g«) > or < 0, we write a« = 1 or = 0, 
respectively. Binary positions for the whole chain can be compactly encoded in a binary 
integer (signature) 



Thus we have defined a partition of 2iV-dimensional phase space cut by N hyper-planes 
qi+i = Qi into 2 N cells labeled by signatures S G {0, 1, ... 2 N — 1}. By a slight abuse 
of language, we will use a term signature S also to refer to a phase space cell denoted 
by S. If the transitions between various signatures are rare, which is obviously the case 
if the potential barrier AU is high, it seems to be meaningful to concentrate on the 
(statistical) dynamics of transitions between the signatures. For the periodic boundary 
conditions the coordinate differences must fulfill the constraint Y^Li (<7i+i — g?) = 
which in the limit of infinitely high potential barrier and even N translates into the 
condition YliLi a i = N/2. For odd N the number of pairs in the left and the right well 
should differ by 1. This presents only unnecessary technical complications and from 
now on we will assume N to be even. If the barrier is finite, signatures with different 
number of pairs in the left and right well are possible to visit. We will call signature S 



N 

S = (a N . . . a 2 a 1 ) 2 = X! a « 2 



(3) 



i=i 
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to be of order i if J2iLi a i = N/2±i. The number of different signatures Mj of order i is 
^ JV \_ N\ 
vN/2-i) " 



Mi 



(4) 



) (N/2-i)\(N/2 + i)V 

For sufficiently high barrier the system will spend most of its time in signatures of order 
0. We will specifically concentrate on signatures of order and treat signatures of higher 
(mostly 1st) order only as 'tunnels' for transitions between different signatures of order 
0. Figure |2| shows an example of a transition between two order signatures via an 
intermediate short lived signature of order 1. First we have to know for which values 
of parameter a and size N will the description by the signatures of order only be 
adequate, i.e. will the relative time spent in higher order signatures be negligible. 
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Figure 2. Time dependence of differences (qi+i(t) — qi(t)) for N = 14, a — 4.8. For 
time t < -20 the system is in the order state S = 01010110101010 2 , at t « -20 
the system jumps to the temporary order 1 metastable state 5" = 010101 IOOOIOIO2 
and shortly thereafter at t = jumps into the new stable order state S' = 
011101 IOOOIOIO2. Such a transition from S to S' will be called a jump of length 
d — 6. Dashed lines denote the position of the saddles qi+i = q%- 



3.1. Estimating the fractional volumes of higher order signatures 

We will now make a rough theoretical estimate for the time spent in signatures of 
different orders. For an ergodic system this time is proportional to the measure (volume) 
of the phase space covered with the signatures in question. Let us designate by t{ the 
time spent in signatures of order i. We have 

t (E) oc [ dS dVi ^ E ) (5) 
J order i \VH\ dE ' w 
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where the region of integration is the part of the energy surface H = E intersecting with 
the full set of signatures of order i and Ti(E) is the total volume of the phase space in 
the signatures of order % and with energy less than E. First, let us make an estimate 
for T (E). Potential around both minima is to the lowest order harmonic, therefore we 
can roughly say Vq(E) ~ T(E) N , where T(E) is the volume for one harmonic oscillator, 
that is T(E) = 2-nEj^fa. For the signatures of order i the argument is very similar. In 
this case there are 2i particles more in one potential well than in the other. Because the 
trivial condition Y,i=i — Qi) = must still be satisfied, the equilibrium positions 
of pairs will not be at ±y/oi any more but will be shifted for Aq in order to keep the 
center of mass of all pairs at 0. This gives for the shift Aq = 2i\fa/N and contributes 
£i = aAq 2 = 4i 2 a 2 /N 2 to the energy density, provided the cubic part of the potential is 
negligible, i.e. 2i/N <C 1. Harmonic approximation for the volume Ti(E) can therefore 
still be used, but now at energy E — Nsi. Since we are interested only in the dependences 
for large N we can omit ('cancel') energy derivatives coming from ([5]), so that we have 

* r(E - e ' ]N <1-«Y -evf^g-). (6) 



to r( £ )» V e 

We have fixed e = 1 so that the relative fraction of signatures of order 1, t\/to, falls 
exponentially in a 2 /N. We must stress that the exponential dependence is expected only 
when 4a 2 /N 2 <C 1 and that the overall prefactor in front of an exponential function 
could still depend on N. This is nicely confirmed by the numerical data in figure []. 
Grouping the signatures by their order is meaningful since members of each group have 
the same phase space volume and the same minimal potential energy. The minimal 
potential energy of order signatures is Eq = 0, of order 1 is Ex — 4a 2 /N 2 and of order 
% is Ei m i 2 E\. The fraction of time spent in higher order signatures is just a power of 
relative time spent in order 1 signatures, (ti/t ) ~ (tx/t o y 2 and this can be reduced by 
increasing a. 

Of course, for the transitions to be possible at all, the barrier height must be smaller 
than the total energy E. This yields the condition 

a 2 1 '4 < N (7) 

which is just the opposite of the condition to keep the tx/t small. Nevertheless, we 
can still reduce the fraction of higher order states to around exp (—16) ~ 10~ 6 , which is 
small. But the condition ([/]) can be problematic for small N. Indeed, for the smallest 
nontrivial case of N = 4 we have numerically found that there are no transitions between 
signatures for a > 2.8. For smaller a we do find transitions, but there the system exhibits 
non-ergodic behaviour. Condition (|7|), although fulfilled, seems to be too weak and is 
still keeping the phase space cut into non-connected parts. The smallest system worth 
studying regarding transitions is therefore N = 6. 

In all numerical work a and N have been chosen so as to keep the fractional volume 
of higher order signatures small. We have therefore studied the dynamics on a set of Mq 
signatures of order exhibiting transitions thru short-lived order 1 signatures as shown 
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Figure 3. Fraction rj — (SSi*<)Ao ~ ti/to of time spent in signatures of higher 
order for AT = 6, 8, 10, 12, 16, 20 and 30, from right to left, respectively. Lines connect 
points with the same value of N. Inset shows by iV-dependent prefactors rescaled 
values of r\ and exponential function (0) (chain straight line). 



in figure |2|. The most important physical scale that characterizes such a transport is the 
average time r between subsequent transitions. 

3.2. Ergodicity and average transition time r between signatures of order 

Though the main body of this paper is concerned with the transport and decay of 
time-correlations, i.e. with the system's mixing property, we should first mention, that 
we have also checked the weaker dynamical property of ergodicity directly. We have 
employed a method of Robnik et al |T| of comparing the rate of visiting of different 



phase space cells (in out case, order signatures) with the rate for a fully random 
dynamics. The results turned to be fully consistent with (uniformly) ergodic behaviour 
of our IFPU model in the high-barrier regime discussed above. 

Now we define the time scale r to be an average time from the point when the orbit 
enters a certain order signature S to the point where the same orbit enters the next 
order signature S' (S' ^ S) concluding the transition from S to S'. The average is 
taken over many different orbits with microcanonically distributed initial conditions, or 
over one very long orbit if we assume ergodicity. We obtain an approximate functional 
form for the dependence of r on N and a by similar arguments as for the relative 
volume of higher order signatures. The probability for a jump will be estimated by the 
ratio between the phase space volume of the set of states just before a jump I\ and 
the volume of the set of equilibrium states r eq . The approximate volume T eq (E) of the 
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phase space of equilibrium states has been derived in the previous subsection and is just 
T eq (E) ps T(E) N = (2tc I \/a) N E N . In the state just before the transitions there should 
be more than AU of energy in one pair of particles. This energy will be denoted by aa 2 , 
where a is for now some unspecified constant. Certainly a must be greater than 1/4, 
but because our oscillators are coupled we allow for the possibility that the energy of a 
pair must actually be bigger than the barrier height. Phase volume T t (E) for a state 
with one pair having the energy aa 2 and the other N — 1 pairs residing in the minima 
is 

r-E-ac?- ( Ovr \ N ~ X 

T t (E) = J dee^l^L) (N-1)(T(E)-T(e)) 




- aa 2 ) 1 "- 1 . (8) 

Probability for a jump is proportional to the quotient of derivatives dT t (E)/dT eq (E). 
Reciprocal average time between transitions l/r(N,a) is proportional to the above 
quotient and to the frequency of oscillation around potential minimum. Using also 
E = N, we finally get 

t(N, a) ps ^^jz , (9) 

V« acr (1 — aa z /N) N 



where A is some overall numerical factor. We have determined numerical parameters a 
and A by fitting the numerical data for different lengths iV in the range N — 6, . . . , 30. 
Both parameters do depend on N for small chains, but are, of course, independent 
of a. The value of A is between 0.02 and 0.05, whether the a has dependence 
a = 1/4 + 0(1/N). In the limit of large chains a has the value 1/4 predicted by 
our simple physical picture. This picture can also be confirmed by looking at the local 
energy density just before the jump. Indeed, local energy of a pair is exactly aa 2 for all 
chain sizes. So, the parameter a is not just some fitting parameter in equation @ but it 
has a clear physical interpretation. In figure |] we depict the estimate of r (H) together 
with the numerical data. Worth noticing is also an interesting scaling law, although only 
approximate. Global dependence of r(N,a) is such that we can write approximately 

t(N, a) PS /(4A ^* 2) exp (« 2 /4), (10) 
'a 



where f(x) is some parameter-free function. The shape of this function can be seen in 
figure 0. 

We have also numerically computed statistical distribution of transition times r n , 
i.e. the time intervals between entering different subsequent signatures of order such 
that r = (r n ), for a fixed N and a. This distribution is shown in figure |6] and is clearly 
seen to be exponential. There is a discrepancy only for extremely small times. This is a 
consequence of the duration of intermediate unstable order 1 state that occurs between 
order states (see figure |2|). The fraction of order 1 states was in this case ti/t = 0.0013 
which is the same as the time scale of the discrepancy in the figure. Our distribution 
should be a convolution of distributions of to and of t\. But the precise explanation of 
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Figure 4. In the top figure experimental values of r for N — 6, 8, 10, 16, 20 and 30 
(symbols from left to right) and semi-theoretical estimate (with continuous curves) for 
r ([)]) are drawn. Errors in numerical r are of the same order as the symbol sizes. In 
the bottom figure we plot the dependence of r on r\ = (Y]'^_ 1 i,)/io- The data and the 
symbols used are the same as in the top figure, while the straight line segments are 
here merely connecting the points with the same value of N. 



this discrepancy is the following: we obtain too many of "short" jumps because of the 
way we measure the time when the system has arrived into a new order signature. We 
mark the system as arriving into a new signature exactly at the top of the barrier, when 
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Figure 5. Numerical data for r(N,a) divided by exp (a 2 /4)/y / a versus AN /a 2 . The 
solid curve is only drawn to guide an eye, perhaps suggesting some universal scaling 
function. 




t/x t/t 

Figure 6. Probability distribution of times t = t„ (in a relative scale t/r) between 
two different consecutive order signatures. On the right we plot enlarged distribution 
for small times. All is for the chain of length N = 20 and a — 5.4 corresponding to 
r = 3500. Dashed line is an exponential function. 



the difference — changes the sign and not when it relaxes to the bottom of the 
well. This allows for some 'fake transitions' in the regions of phase space where two or 
more order signatures almost touch with a saddle-like region in between (for instance 
near coordinate origin = 0). There it is possible that the system, which although it 
gets over the barrier, subsequently relaxes into some nearby well because it does not 
have the 'right' momentum vector. 
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4. Markov process and its transition matrix 

We decided to describe dynamics with a sequence of binary signatures instead of a full 
orbit x(t) with the hope of some simplification. We are hence looking for some simple 
statistical description of the transitions between signatures. Stochastic system is fully 
described by giving conditional transition probabilities P(S, S', S", ...;£,£',.••)■ The 
simplest kind of a stochastic process is called Markov process. For a Markov process the 
probabilities P(S, S'; t) contain entire information about the system. We will shortly 
write a matrix P(t) instead of P(S,S';t). For now let us define P(S,S';t) on the full 
space of 2^ signatures although we will later be interested only in transition probabilities 
among M signatures of order 0. The transition matrix P must satisfy the condition 

p(t + t') = p(*)P(0- (ii) 

The matrix P(i) contains therefore some redundant information. We can write instead 

P(t) = exp(Ft), (12) 

where for a Markov process the matrix F is time independent. The matrix element 
F(S, S') is a probability flux from the signature S to the signature S", and can be 
defined by the limit 

F = lim P(t) ~ I . (13) 
t^o t 

The conservation of probability imposes the condition J2s' F(S, S') = for the matrix 
elements. 

By a statistical description, using P(t) instead of x(t), we lose information about 
the details within a single cell of a partition. In other words, we study the dynamics 
on the space of functions that are constant over one phase space cell. This can be 
formalized by introducing characteristic functions P>s(x) over signatures S 

B s {x) = { 1 X f S . (14) 
v ; \ otherwise v ; 

Characteristic functions Bs span the space of all observables which are piece-wise 
constant over the signatures. Every piece-wise constant observable W(x) can be 
expanded over the base functions B s as W(x) = J2s w sBs( x ) with the expansion 
coefficients ws given by 

where the brackets ()e denote a microcanonical phase space average over the energy 
surface H = E. The matrix elements of a transition matrix P are nothing but the 
correlation functions between the characteristic functions 

P(S,S';t) = (B s ,(t)B s ) E /(B s ) E . (16) 

As the transition matrix propagates the probabilities between the signatures we can 
write the time dependent vector w(t) = (ws(t),ws'(t), . . .) as 

w(t) = P(t)w(0) = exp (Ft)io(O). (17) 
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Similarly, the autocorrelation function of the observable W is 

(W(*)W(0)) E = w T exp (Ft)w. (18) 

We can see that the behaviour of correlation functions is determined by the spectrum 
of the (Markov) propagator exp (Ft) . 

From now on we assume that the fraction of higher order signatures is very small 
r) ~ ti/to <S 1 and we restrict the matrices P(t) and F on the Mo-dimensional subspace 
of signatures of order 0, which essentially contain all non- vanishing matrix elements. 
Numerically we calculated the transition matrix P(t) from one very long orbit as an 
average 

P(5,5';t) = ^^%^, (19) 

\0S,S(t'))t' 

where the brackets () t / denote a time average over the orbit and 6s,s(t) has value 1 if 
the orbit is in a signature S at time t and otherwise. For small t the average can be 
performed, and by using expression ([Uj), yields for F(S, S'), (S ^ S') 

F(S,S>) = ^±, (20) 

ts 

where n(S, S') is the number of direct transitions between the order signatures S and 
S' and ts is the total time spent in the signature S. Here a small comment regarding 
higher order signatures is in order. We saw in figure ^] that between each order 
states there is a short intermediate order 1 state of duration ~ rt\/to- So there are no 
direct transitions between order states and all fluxes F(S, S') (|20f) among M order 
signatures are strictly zero. As we have decided to study only (indirect) transitions 
between order states (for the parameter values where the fraction of higher order states 
is negligible) we must somehow circumvent this difficulty One solution is to calculate 
the derivative of the transition matrix P(t) in equation (|T3|) not at time t — but at the 
time t with rti/t <^ t <^ r. In this case the number n{S, S') in ( p0|) is the number of 
jumps between order states S and S' with one short intermediate order 1 state. Still 
more simple and convenient solution is to replace the orbit x{t) by a sequence of times 
ti when the orbit hits the boundary of a cell of order coming from outside. Orbit x(t) 
is therefore replaced by a sequence of pairs . . . , (Si,ti), (Si + i,ti + \), . . . which tell you 
that at time ti orbit x(t) came to the order signature Si, then until time ti+\ was in 
signature Si or any higher order signature and at time t i+ i arrived into the next order 
signature S i+ i (S i+ i ^ Si), etc. Then we define n(S, S') as the number of subsequent 
pairs (S,S') in the sequence {Si}, and ts = J2i 1 ~ S {U+i — U), hence the flux matrix is 
calculated as 

F(S, S') = r- (21) 

This definition is used in the numerical calculation of this paper. 
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4-1. Tests for the Markov process 

Before describing the process of transitions as Markovian we must check if this is at 
all permissible. Checking the rigorous conditions for the partition to signatures to 



be Markovian, see e.g. ||12|| , is a considerable mathematical problem, at least to our 
judgement. On the other hand, heuristic arguments and numerical results support the 
idea, that for sufficiently large transition times r the process will indeed be Markovian. 
The system has positive Lyapunov exponents.|] Vaguely this means that the system 
quickly "forgets" its history. If the average time between the jumps r is longer than 
the Lyapunov time, we can expect transitions to be Markovian. However, the most 
straightforward test is to check explicitly the composition formula (fLl|). 

Condition P(2t) = P 2 (t) is trivially fulfilled for small and for large times t, t <C r 
and t>r, respectively. We checked it for time t = r by calculating the quantity 

||P(2t)-P 2 (t)|| e 
= l|P(2r)|| E ■ (22) 

where ||A||e = (Y^i.j \^i,j\ 2 Y^ 2 is the Euclidean norm of a matrix A (all results have 
also been re-checked by using the spectral norm giving almost identical results). As the 
matrix P(i) will be calculated numerically from one very long but finite orbit it will 
have statistical error ctt due to the finite number of simulated transitions. In addition, 
a will also have a contribution o"m from a systematic error because the process may not 
be precisely Markovian. We have assumed a 2 = + o~\ and made an estimate for 
the statistical error ctt as a square root of a number of average transitions per matrix 
element n. It is 



1 N M ° 

aT = c " c TV^o' (23) 



where Nq is the total number of all transitions between the signatures of order and c 
is some unspecified numerical constant which has been determined from the numerical 
data. Its value turned out to be c « 0.7. Results for the dependence of a systematic 
error <tm are shown in figure [7[ The error <7m(t) goes to zero with increasing r, as 
expected. Even more interesting is the fact that the error also seems to decrease with 
increasing N at constant r. This is interesting as it means that with increasing N we 
do not have to keep the fraction of higher order states small in order for the process to 
stay Markovian. @[] 

Whether a given system can be described as a Markov process can be checked also 
by calculating the transition matrix P(t) and comparing it with the exponential formula 
([Hp as a function of time. P(t) can be calculated numerically directly from definition, 
for instance by using the formula ([19]). This definition would have a meaning also if the 
transitions were not Markovian. We then compare the correlation functions between the 

I We have numerically checked that the distribution of Lyapunov exponents is approximately linear, 
as is typically the case for sufficiently chaotic systems |l3| . 

% In such a case one could bring higher order signatures back and describe a complete hierarchy of 
transitions between signatures of various orders as a Markov process on 2 N dimensional space. 
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Figure 7. Dependence of cm on r in figure (a) and on 77 = (X^i U)/to in figure (b). 
Pluses, crosses, stars and squares are for the N = 6,8,10 and 12, respectively. 



characteristic functions ( [TBI) and the matrix elements of the propagator (S\ exp (Ft)|S"). 
These calculations should agree only if the flux matrix F(t) = logP(t)/t were time 
independent and this would signal that the transitions are described by a Markov 
process. Numerical results for this test are shown in figure |8|. Correlation functions 
calculated in both ways have been computed for the characteristic function B$ = \S) on 
the signature S = 54 and for the macroscopic observable with the vector ws = S mod 2, 
i.e. characteristic function on a union of half of all signatures with the first bit a± = 1, 
denoted by [si). It can be seen from the figure that the autocorrelation function 
of the observable | si) begins to fall as exp (—Ait), where Ai is the biggest nontrivial 
eigenvalue of F, very early. This is a consequence of relatively large overlap between 
the macroscopic state | si) and the eigenstate |i>l) corresponding to the eigenvalue 
Ai. We will show in subsection ^ that |(slK)| 2 ~ 1/N while |(SH)| 2 ~ 1/M , 
so |(slh)| 2 > |(%i)| 2 (for iV> 1). 

5. The form and the spectral properties of the flux matrix 

Something can be said immediately about the spectrum of F. Conservation of 
probability implies all eigenvalues Aj to be smaller or equal to zero. Further, because the 
system is believed to be ergodic, there should be a non-degenerate constant eigenvector 
with the corresponding eigenvalue Ao = 0. This is all the consequence of the Perron- 
Frobenius form of the matrix P. For finite N the spectrum of the flux matrix is discrete 
and the correlation functions will fall asymptotically as exp (—Ait) where Ai is the largest 
nonzero eigenvalue. An interesting question is what happens in the thermodynamic limit 
iV — > 00. If there is a spectral gap we will have an exponential decay, otherwise some 
anomalous behaviour can occur. 
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Figure 8. Comparison of the correlation functions calculated from the propagator 
exp(Ft) (dashed curves) and from directly determined transition matrix P(S,S';t) 
(full curves). In figure (a) are autocorrelation functions S = S' and in figure (b) 
cross-correlation functions S ^ S'. Bottom pair of curves in figure (a) is for the 
signature S — 54, the middle pair for the one particle state \sl) and the top curve 
is a simple exponential decay with the biggest nonzero eigenvalue of the matrix 
F. From all correlation functions we have subtracted equilibrium value which is 
P(S,S;t — > oo ) = 1/Mo and 1/2 for the state \sl). In figure (b) are correlation 
functions between the signature S — 54 and the signature S' = 58, S' = 116, and 
S' = 135, from top to bottom, respectively. The chain length and the parameter were 
N = 8, a = 3.83, r - 1000. 



5.1. Numerical calculation of the flux matrix 

Non- vanishing matrix elements of F are expected only between the signatures that differ 
just by two bits provided that ti/ta <^ 1. In other words, two bits Oj = 1 and aj = only 
switch their position. Such a transition will be called a (bit) jump of length d, where 
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0.80 


0.74 
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0.78 
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0.82 


0.76 


0.82 


0.86 



Table 1. Dependence of Cd on distance d for different r and iV. 



^ = |j — ^1 is the distance between the two bits involved. Example of a such transition 
has been shown in figure We have numerically calculated complete matrices F for 
sizes N = 6, 8, 10 and 12 and different r (or a), in order to check the above hypothesis 
and to obtain some knowledge about the sizes and the structure of matrix elements. 
The only matrix elements which are really non-vanishing are those of the bit jumps, 
and what is more, the value of the matrix element depends only on the distance d of 
the jump and not on the particular signatures involved. This practically means that in 
addition to the diagonal elements we have only N/2 different matrix elements in the 
flux matrix F. This number must be compared to the total number of matrix elements 
Mq which grows exponentially with N. We have checked this also for the matrices F of 
sizes up to N = 30 but in this case we have compared only the jumps of length d with 
different number of bits 1 on the sites between the two jump sites. We confirmed that 
the probability fluxes of jumps depend only on the length d. It is therefore meaningful 
to introduce dimensionless coefficients q which are ratios between the matrix elements 
of a jump of length d and a jump of length d = 1. Therefore c% = 1, by definition, and 
the rest of N/2 — 1 coefficients q together with the average transition time r is all that 
we need in order to specify the flux matrix F completely. We can write 

F = I c, (24) 

where the matrix C is just a matrix involving the coefficients q only. Numerical values 
of Cd for different N and r are listed in table [I]. It can be seen that for not too big r 
(e.g. r = 11 and N = 10) the coefficients q decrease monotonically with the distance 
d. For bigger r (e.g. r = 621 and N = 10), this dependence ceases to be monotonic 
but becomes slightly well-shaped with a minimum at around d ~ N/4. There also 
seems to be a general trend that with increasing N, at a constant fraction 77 pa t±/to, the 
coefficients Cd all approach 1. This can be explained by the following heuristic argument. 
To keep the fraction ti/t of higher order states constant with increasing iV we must 
increase the barrier height AU. As a consequence, the transition time r increases and 
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correspondingly also the average time rti/t spent in the non-equilibrium intermediate 
state of order 1 in between the states of order 0. At a constant ti/t the transition 
time r depends approximately exponentially on the chain length N (see figure |4j). On 
the other hand, the relaxation time for the non-equilibrium energy distribution at the 
moment of transition grows perhaps only linearly with iV (inversely proportional to 
the sound speed), but certainly slower than exponentially. For sufficiently large N the 
relaxation time will therefore be smaller than the time spent in the intermediate order 1 
cell. Initial locally non-equilibrium distribution at the beginning of a jump will therefore 
relax before the transition to the new order state, which will therefore all be equally 
probable. Because of this we can assume that the coefficients q will not depend on d in 
the thermodynamic limit. This hypothesis still needs further verification as the limit of 
high N could not be tested due to the fast growth of the dimension Mo with N. Despite 
of that, the numerical results are consistent and point in the right direction (see table 
|]). We will show now, that our flux matrix has a specially appealing form which can 
be described using a formal connection to a certain Hamiltonian of a quantum spin-1/2 
chain. Based on this we will also deduce the essential spectral properties of the flux 
matrix. 



5.2. Correspondence between the Markov process and the quantum spin chain 

We can make the announced correspondence by connecting the eigenvalue problem for 
the flux matrix F with the eigenvalue problem for the quantum Hamilton operator H over 
Mo-dimensional Hilbert space with the same matrix elements as F. This correspondence 
would generally be of no particular use since the Hamiltonian H does not represent any 
simple quantum system. In our example of IFPU chain this is not the case, since 
F has a particularly nice and simple form, with non-vanishing elements only between 
the signatures connected by a bit jump where signatures are sequences of binary bits, 
and 1. This should immediately remind us of Heisenberg chains of quantum spin- 
1/2 particles. Let us write the Hamilton function for one-dimensional ferromagnetic 
Heisenberg spin chain of length N, with Pauli variables o,j = (<rj, <rj, af) and coupling 
constants J(d) 

1 N 1 N ( 1 



l JV l JV r 

H = "^ E J{d)a r a 3+d = -- £ J(d)\c 
Z hd=l z j,d=l K 



tt+d + o K +(7 7+<* + °i a t+d) \, (25) 



where periodic boundary condition is assumed cr N+ j = <jj and J(d) = J(N—d), J(0) = 0. 
Operators af = aj ± iaj are the standard raising and lowering operators. Quantum 
state of the spin chain, an eigenstate of cr|, will be denoted by the signature \S), with 
the obvious interpretation. Bit aj = 1 or denotes spin j up or down, respectively. 
Now we can calculate the matrix elements of the Hamilton operator (|25|) . For the off 
diagonal elements we get 

f —2J(d) if there is a jump of length d connecting S and S' 
(S\H\S) = < , ta . (26) 

x 1 1 7 otherwise v ; 
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Diagonal elements are 

N/2-1 

(S\E\S) = -J(N/2) [N/2 - 2s N/2 (S)] ]T J(d) [N - 2s d (S)}, (27) 

d=l 

where Sd{S) is a number of different signatures that can be reached from the signature 
S with a jump of length d. For any signature S we can write an identity 

N / 2 , N 



d=l 

T2 



E **(S) = ( y ) • (28) 



By denoting s = MqN / (2N — 2), the following equality is also valid 
^ f s d=l,... N/2-1 

^ Sd{Sl) = { s/2 d = N/2 ■ (29) 

The ground state for the ferromagnetic Heisenberg Hamiltonian (|25|) can be found 
immediately for any J(d) > 0. We have to keep in mind that the Hilbert space is in our 
case spanned just by M signatures of order and not by all possible signature states 
as is usual for the quantum spin chains. Order i of the signature is simply an eigenvalue 
of the S z component of the total spin, namely S z = SjCr|/2 = i. But the operator S z 
commutes with the Hamilton function and the Hilbert space is therefore a direct sum 
of subspaces labelled by eigenvalues of S z . We choose M dimensional Hilbert subspace 
with S z (= %) = 0. Ground state denoted by |0) is 

i m 

|0> = -7jFp£l<Si>> (30) 



with the energy 



'JV/2-1 1 



E = (0|H|0) = —N g J(d) + -J(N/2) J . (31) 

Now that we have the matrix elements of H and the ground state, we can see that 
if we prescribe J(d) to equal Cd multiplied by the constant prefactor in front of the 
matrix C in equation for F fl24|) , we can formally write our flux matrix in terms of the 
Heisenberg spin Hamiltonian (|25|) 

F = -H + E I. (32) 

Transition probabilities P(S, S'; t) can now be written as 

P(t) = exp (Ft) = A(t) exp (-/?H), (33) 

where we wrote (3 = t and A(t) = exp (E t). Expression for P has the same form as 
the density operator for the quantum canonical distribution at temperature 1/t. Time 
dependence of P is therefore the same as the dependence of the canonical distribution on 
cooling. Understanding time dynamics of the IFPU model is equivalent to the cooling 
of the ferromagnetic Heisenberg spin chain or its imaginary time dynamics. In the 
limit t — ► oo the matrix elements of P go towards the constant value 1/Mq and in 
the corresponding Heisenberg spin chain any non-equilibrium distribution relaxes to the 
ground state |0). 
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For higher lying states we numerically solved the eigenvalue problem for the flux 
matrix at various r and sizes up to iV = 12. One such example of a numerical spectrum 
is shown in figure [| We have said that with increasing N, at constant tx/to, we 
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Figure 9. Spectrum of F for N = 8 at three different r = 17, 180 and 1030 with 
stars, pluses and crosses, respectively. With boxes is plotted the spectrum in the case 
c<j = 1 ( |35| ) for comparison. 



expect coefficients q to approach a constant value q = 1 independent of d. This 
case corresponds to the simplest Heisenberg spin chain with uniform coupling. Apart 
from the constant factor the Hamilton function ( p5[ ) is in this case H = — 2S 2 + 3N/2 
and the matrix F reads 

2 ( _^ iV(iV + 2).' 



2S Z 



I 



(34) 



N 2 t \ ~ 2 

where S 2 = 1/4(1?! + ... + a^) 2 - Eigenvalues of S 2 are S(S + 1), where 5* is the 
quantum number of the total spin with values from S — 0, . . . , N/2. The Hilbert space 
is composed of order states with S z = 0. If we denote by Xj the eigenvalues of operator 
(|34|) and by n(j) the corresponding multiplices we have 

Aj(N+l-j) 



A, 



n{j) 



(N 



N 2 r 
1 - 2j)N\ 



j = 0,...,N/2. 



(35) 



(N + l-j)\(j)\ 

Particularly interesting are the first two eigenvalues Ao and Ai. The biggest, non 
degenerate eigenvalue Ao = belongs to the ground state |0). The largest nontrivial 
eigenvalue is Ai = —A/Nt with the multiplicity n(l) = N — 1. In the thermodynamic 
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limit and keeping r constant, it goes to zero as 1/N. The correlation functions therefore 
decay slower than exponential in the thermodynamic limit. From the figures ^ and 
10| it can also be deduced that with increasing r the spectrum is indeed approaching 
the spectrum for the q = 1 (pop, as predicted. What is more, the eigenvalue Ai 
seems to be monotonically approaching the limiting case (c<j = 1) from above. We can 
therefore conclude by observation that the correlation functions of observables which 
can be spanned by \S) decay slower than exp (— 4t/iW). 
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Figure 10. Enlarged first multiplet for N = 1Q and r = 11,40,250 and 620 with 
pluses, squares, full squares and crosses, respectively. Referential degenerate set of 
first N — 1 eigenvalues at c<j = 1 is plotted with stars 



Now that we expect the spectrum to be close to the one for Cd = 1 we can 
explain the decay of the correlation functions for one particle state | si) and the one 
signature state \S) in figure j|. State vector for \S) has only one component different 
from 0. Eigenvector \vl) corresponding to the eigenvalue Ai has on the other hand 
all components approximately of the same order, i.e. 1/a/M . Square of the scalar 
product is therefore ((S^fi)! 2 ~ 1/Mq. One particle state | si) is proportional to the 
\sl) ~ erf erf 1 0). Because we act on the state |0) with the eigenvalue S = N/2 with the 
operator er^af , we can decompose the state | si) as a linear combination of states with 
the eigenvalue S = N/2 and S = N/2 — 1. State | si) can therefore be written as a sum 
of iV states, one of which is also \vl). If we assume all the expansion coefficients to be 
of the same order we immediately obtain |(sl|t>i)| 2 ~ 1/N. 
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6. Discussion and conclusion 

We divide the phase space of a iV-body Hamiltonian, namely the inverted FPU model, 
into 2 N cells that are uniquely tagged by a binary number, called a signature (Q). 
Signatures are ordered according to the absolute difference in the number of particle 
pairs/bits in the left/0 and right/1 potential well. For a sufficiently high potential 
barrier an approximate Markovian description on space spanned by order signatures 
is possible. The accuracy of the Markovian property has been checked numerically to 
be increasing with increasing chain length N at a constant fraction of higher order 
states. In the flux matrix F we have non-vanishing matrix elements only between the 
signatures connected by an exchange of a pair of bits. This enables for the formal 
correspondence between the IFPU and ferromagnetic Heisenberg quantum spin-1/2 
chains. Understanding time dependence of a Markov transition matrix is equivalent to 
cooling or imaginary time dynamics of a quantum spin chain. If we increase the chain 
length N at a fixed fraction of higher order states, the transition rates are expected 
to become independent of a jump length d. Such trend is confirmed by a numerically 
calculated matrices F. In this limit exact eigenvalues of a flux matrix can be obtained 
explicitly. The largest nontrivial eigenvalue in this case is Ai = — 4/iVY and goes to in 
the thermodynamic limit. There is also a numerical evidence that the largest nontrivial 
eigenvalue in a general case (for a non-constant transition rate) is strictly larger than 
— 4/iVY. Therefore, the correlation functions of the observables which can be spanned 
by functions \S) decay asymptotically slower than exp (—4t/Nr). 

Unfortunately these results do not imply directly the behaviour of more general 
correlation functions, such as current-current correlation which is needed in order to 
understand the transport properties. We have made some numerical calculations of 
current- current time correlation functions of IFPU chain for different sizes N and 
different parameters. As expected, the decay of correlations for observables that vary on 
time scales smaller than r (e.g. current) is faster than the decay of correlation functions 
spanned by the piece- wise constant basis \S). Transition from the algebraic (anomalous 
conductivity) to the exponential (normal conductivity) decay of current autocorrelation 
function occurs rather abruptly at a ~ 3 (for iV > 20). 

Finally, though we suggested several interesting properties of the above model, we 
must admit that most of our conclusions are based on numerical evidence. Therefore, 
we believe that it should be a challenging and not impossible future task to try to 
provide more rigorous justifications (and perhaps proofs) of our results. Establishing 
rigorous asymptotic Markovian property would enable one to systematically code and 
enumerate all the many-body (unstable, hyperbolic) periodic orbits and to use them 
explicitly in a classical or semi-classical trace formulae, for example to calculate the 
transport coefficients directly. We feel that IFPU chain may become a useful toy model 
of a chaotic field theory |14| . 
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